Introducción

El análisis espacial y espacio-temporal es una extensión natural a las capacidades de R. Existen múltiples paquetes para manipular y analizar datos espacio-temporales en R. En esta clase vamos a dar una introducción a esto.

matrices y data frames

Recordando, una matriz es un arreglo bidimensional con una cierta cantidad de renglones y de columnas, por ejemplo un cuadrado de \(3*3\) que contenga los números del 1 al 9

matriz <- matrix(seq(1,9),nrow=3,ncol=3)

matriz
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

Un data.frame es también un arreglo con la diferencia de que puede contener más de un tipo de datos. Esto es, puede alojar texto, números y valores booleanos simultáneamente.

numeros <- c(1,2,3)

texto <- c("hola","como","estas")

booleanos <- c(TRUE,FALSE,TRUE)

data_frame <- data.frame(numeros,texto,booleanos)

data_frame
##   numeros texto booleanos
## 1       1  hola      TRUE
## 2       2  como     FALSE
## 3       3 estas      TRUE

Ambos tienen asociadas coordenadas de la forma \((renglón,columna)\) que comienzan de arriba a abajo y de izquierda a derecha. Estas permiten recuperar los valores que guardan, por ejemplo la entrada \((renglón,columna)=(2,1)\) para la matriz y el data.frame que definimos antes corresponden a:

matriz[2,1]
## [1] 2
data_frame[2,1]
## [1] 2

También se pueden hacer búsquedas condicionales sobre estos objetos, esto es escribir en código peticiones como: dime qué valores de la primera columna de la matriz son mayores a 1.

Cuando se deja vacía una dimensión, R entiende que le estás pidiendo que incluya todos los valores existentes en ella. Por tanto, la pregunta ¿qué valores de la columna 1 de la matriz son mayores a 1? se expresa de la siguiente manera:

condicion <- matriz[,1]>1

condicion
## [1] FALSE  TRUE  TRUE

La primer columna está compuesta por los valores 1 (que claramente no es mayor a 1), 2 y 3 por lo tanto el resultado es FALSE, TRUE, TRUE.

La misma manera de hacerle “preguntas” a una matriz se puede lograr con un data.frame

condicion <- data_frame[,1]>1

condicion
## [1] FALSE  TRUE  TRUE

Algunas preguntas, por la naturaleza de los datos no tienen sentido. Por ejemplo de las palabras “hola”, “como”, “estas” ¿cuáles son mayores a 1?

Rasters, primeros pasos

Un raster es simplemente una matriz de números qué tiene la intención de usarse como base para generar imágenes. A estos se les puede asociar una estructura espacial.

En R, existe el paquete “raster” que está ampliamente desarrollado y se recomienda como espina dorsal para cualquier análisis que incorpore imágenes de satélite.

La mayoría de las imágenes de satélite se pueden cargar usando la función raster(). Esta misma permite insertar una matriz en un objeto raster de R.

library("raster")
## Loading required package: sp
## Warning: no function found corresponding to methods exports from 'raster'
## for: 'overlay'
matriz_raster <- raster(matriz)

matriz_raster
## class       : RasterLayer 
## dimensions  : 3, 3, 9  (nrow, ncol, ncell)
## resolution  : 0.3333333, 0.3333333  (x, y)
## extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 1, 9  (min, max)

Como podrá notarse, el objeto matriz_raster es ya un raster. Si bien, uno poco interesante y sin proyección.

Se mencionó que los rasters están pensados para ser la base de imágenes. El paquete “raster” contiene funcionalidades básicas de visualización de objetos raster.

plot(matriz_raster)

Es muy importante observar que un objeto raster tiene asociados dos conjuntos de coordenadas. El objeto matriz_raster anterior se generó a partir de una matriz de dimensión \(3*3\) por lo que está conformado por 9 píxeles que se ubican a través de su (renglón,columna). Por otro lado, define un espacio abstracto continuo localizado en el extent: 0, 1, 0, 1 (xmin, xmax, ymin, ymax). Por tanto la resolución (tamaño de cada píxel) del raster es de \(0.333*0.333\).

En este contexto, los puntos son objetos espaciales sin área. Por lo que uno se puede localizar en cualquier lugar del espacio que definimos, por ejemplo en las coordenadas \((0.25,0.75)\)

plot(matriz_raster)

points(0.25,0.75,pch=21,bg="red",cex=2)

Utilizando el paquete “rasterVis” que es uno cuya intención primordial es visualizar rasters y sus análisis podemos visualizar este raster muy simple con los valores correspondientes a cada pixel. Esto permitirá observar el resultado de los siguientes procesos que llevaremos a cabo.

library("ggplot2")
library("rasterVis")
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## 
## Attaching package: 'latticeExtra'
## 
## The following object is masked from 'package:ggplot2':
## 
##     layer
gplot(matriz_raster) + geom_tile(aes(fill=values(matriz_raster))) + 
                  scale_colour_brewer(palette="Blues")  +
                  coord_equal() + 
                  geom_text(aes(label=sprintf("%02.0f",values(matriz_raster))),color="white",size=5)

Análogamente a como sucedía con la matriz, podemos recuperar el valor del raster en la entrada \((renglón,columna)=(2,1)\)

matriz_raster[2,1]
##   
## 2

Un raster también es entendido en R como un vector de arriba a abajo y de izquierda a derecha. En este caso el arreglo tiene longitud 9 y la entrada \((renglón,columna)=(2,1)\) es igual a la entrada \([4]\).

matriz_raster[4]
##   
## 2

Lo anterior en combinación con expresiones condicionales nos permite hacer búsquedas y manipular rasters.

Podemos multiplicar el raster en su totalidad por 2 y restarle la constante 1.

matriz_raster_por2_menos1 <- matriz_raster * 2 - 1

gplot(matriz_raster_por2_menos1) + geom_tile(aes(fill=values(matriz_raster_por2_menos1))) + 
                  scale_colour_brewer(palette="Blues")  +
                  coord_equal() + 
                  geom_text(aes(label=sprintf("%02.0f",values(matriz_raster_por2_menos1))),color="white",size=5)

Podemos reemplazar todos los píxeles del raster que cumplan la condición de ser mayores a 4 por 999. Esto se logra expresando la condición de manera análoga a como lo hicimos para matrices. Como un raster se puede ver como un vector, basta con meter la condición en una única dimensión

# reemplazar matriz_raster donde matriz_raster mayor a 4 por 999
 matriz_raster[matriz_raster>4] <- 999

gplot(matriz_raster) + geom_tile(aes(fill=values(matriz_raster))) + 
                  scale_colour_brewer(palette="Blues")  +
                  coord_equal() + 
                  geom_text(aes(label=sprintf("%02.0f",values(matriz_raster))),color="white",size=5)

Hay funciones básicas en el paquete raster que son útiles en mútiples situaciones:

cellStats() calcula un estadístico definido por el usuario sobre un raster de entrada, por ejemplo se puede usar para obtener el promedio de matriz_raster

cellStats(matriz_raster,stat="mean")
## [1] 556.1111

Análisis de Rasters usando R

Como se indicó se puede cargar prácticamente cualquier raster utilizando la función raster(). Muchas imágenes satelitales son en realidad multiespectrales, esto quiere decir que en vez de ser una sola matriz son varias apiladas. Por ejemplo las imágenes RapidEye son imágenes constan de 5 bandas (green,blue,red, red edge, near infrared). Estas imágenes en particular se deben cargar con la función brick() o stack(). Carguemos un recorte de imagen RapidEye localizada en el estado de Chiapas. Si se llama al objeto se desplegarán los metadatos del mismo:

chiapas1 <- brick("/home/jequihua/Documents/curso_r_conabio/1crop.tif")
## Loading required namespace: rgdal
chiapas1
## class       : RasterBrick 
## dimensions  : 964, 1056, 1017984, 5  (nrow, ncol, ncell, nlayers)
## resolution  : 5, 5  (x, y)
## extent      : 722995, 728275, 1789575, 1794395  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/jequihua/Documents/curso_r_conabio/1crop.tif 
## names       : X1crop.1, X1crop.2, X1crop.3, X1crop.4, X1crop.5 
## min values  :     3706,     2469,     1208,      974,      722 
## max values  :    13994,    13601,    12661,    11159,    15283
# visualizar las bandas RGB
plotRGB(chiapas1,r=3,g=2,b=1)

Utilizando la función subset() podemos obtener una o más bandas de una imagen multiespectral. Por ejemplo podemos extraer las bandas red (VIS) y near infrared (NIR).

VIS <- subset(chiapas1,subset=3)

NIR <- subset(chiapas1,subset=5)

par(mfrow=c(1,2))

plot(VIS,main="VIS")
plot(NIR,main="NIR")

Utilizando los mecanismos de algebra de mapas descritos en la sección anterior podemos trivialmente calcular el NDVI de esta imagen.

\[ \begin{aligned} NDVI := \frac{NIR-VIS}{NIR+VIS} \end{aligned} \]

chiapas1_ndvi <- (NIR-VIS)/(NIR+VIS)

plot(chiapas1_ndvi,main="NDVI")

Utilizando algebra de mapas y la función cellStats() podemos estandarizar esta imagen. Le restaremos su media y dividiremos esto sobre la desviación estándar de la misma.

chiapas1_ndvi_st <- (chiapas1_ndvi-cellStats(chiapas1_ndvi,stat="mean"))/(cellStats(chiapas1_ndvi,stat="sd"))

plot(chiapas1_ndvi_st,main="NDVI estandarizado")

Con base en el procedimiento anterior, tratemos de generar un ejercicio simple de detección de cambios. Tomemos una imagen RapidEye del mismo lugar pero tomada un año después que la que acabamos de trabajar. Llevemos a cabo el mismo procedimiento para obtener su NDVI estandarizado.

chiapas2 <- brick("/home/jequihua/Documents/curso_r_conabio/2crop.tif")

VIS <- subset(chiapas2,subset=3)

NIR <- subset(chiapas2,subset=5)

chiapas2_ndvi <- (NIR-VIS)/(NIR+VIS)

chiapas2_ndvi_st <- (chiapas2_ndvi-cellStats(chiapas2_ndvi,stat="mean"))/(cellStats(chiapas2_ndvi,stat="sd"))

plot(chiapas2_ndvi_st,main="NDVI estandarizado, 1 año después")

Ahora generaremos una imagen de diferencias, D, a partir de estas dos de NDVI.

\[ \begin{aligned} D = NDVI_1 - NDVI_2 \end{aligned} \]

D <- chiapas1_ndvi_st - chiapas2_ndvi_st

plot(D,main="Diferencias en NDVI")

Aquí las diferencias positivas significan que el NDVI de la fecha 1 fue mayor al de la fecha 2. Diferencias negativas significan que el NDVI de la fecha 1 fue menor al de la fecha 2. Lo que dificulta el estudio de los cambios es que existe un continuo de magnitudes de diferencias. Naturalmente suponemos que diferencias de magnitud cercanas a \(0\) deben considerarse no-cambios pero no es claro de qué magnitud debe ser una diferencia para poder tomarse como un cambio. Vamos a visualizar el histograma de las diferencias para darnos una idea de la distribución de estas.

hist(D,breaks=100,main="Diferencias en NDVI",xlab="Diferencias",ylab="Conteos")

Parecería que la distribución de diferencias es aproximadamente normal. Se sabe que en una distribución normal, 95% de los valores se encuentran a 2 desviaciones estándar de su media. Tomando esto en cuenta definiremos como una diferencia significativa (cambio) a toda aquella que se encuentre a dos desviaciones estándar de la diferencia media. Podemos utilizar lo aprendido hasta ahora para condicionalmente reemplazar en la imagen de diferencias todo aquella que queremos descartar por ser una no-significativa.

umbral_positivo <- cellStats(D,stat="mean")+ 2*cellStats(D,stat="sd")

umbral_negativo <- cellStats(D,stat="mean")- 2*cellStats(D,stat="sd")

# Recordatorio: se pueden combinar múltiples condiciones utilizando operadores lógicos
#
# El operador | se usa para expresar "o"
#
# El operador & se usa para expresar "y"

D[(D>umbral_negativo) & (D<umbral_positivo)] <-NA

Este es el resultado final de este proceso. En R, es tan secillo escribir un raster como leerlo. Para este propósito se utiliza la función writeRaster()

writeRaster(D, filename="/home/jequihua/Documents/curso_r_conabio/diferencias_significativas.tif", format="GTiff", overwrite=TRUE)
## class       : RasterLayer 
## dimensions  : 964, 1056, 1017984  (nrow, ncol, ncell)
## resolution  : 5, 5  (x, y)
## extent      : 722995, 728275, 1789575, 1794395  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/jequihua/Documents/curso_r_conabio/diferencias_significativas.tif 
## names       : diferencias_significativas 
## values      : -2.487393, 3.558255  (min, max)

Resultado final.

Capas vectoriales, primeros pasos

Una Capa vectorial es una estructura espacial conformada por unidades (puntos, líneas o polígonos), asociada a una tabla de datos.

La librería GDAL es una librería Open Source para leer y escribir formatos raster y vectoriales. El paquete rgdal es una interfaz a esta librería que permite la fácil lectura de Shape files en R.

library("rgdal")
## rgdal: version: 0.9-2, (SVN revision 526)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: /usr/share/gdal/1.10
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
puntos <- readOGR(dsn="/home/jequihua/Documents/curso_r_conabio/puntos.shp",layer="puntos")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/jequihua/Documents/curso_r_conabio/puntos.shp", layer: "puntos"
## with 11 features
## It has 1 fields
lineas <- readOGR(dsn="/home/jequihua/Documents/curso_r_conabio/lineas.shp",layer="lineas")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/jequihua/Documents/curso_r_conabio/lineas.shp", layer: "lineas"
## with 3 features
## It has 1 fields
poligonos <- readOGR(dsn="/home/jequihua/Documents/curso_r_conabio/poligonos.shp",layer="poligonos")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/jequihua/Documents/curso_r_conabio/poligonos.shp", layer: "poligonos"
## with 3 features
## It has 1 fields
par(mfrow=c(1,3))

plot(puntos)
plot(lineas)
plot(poligonos)

Cargado un shape file en el espacio de trabajo de R, este se convierte en un data.frame espacial. Por lo que la mayoría de las funcionalidades disponibles para un data.frame aplican.

# tipo de objeto
class(puntos)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
class(lineas)
## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"
class(poligonos)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
# cabeza de los datos
head(puntos)
##   id
## 1 11
## 2 10
## 3  9
## 4  8
## 5  7
## 6  6
head(lineas)
##   id
## 0  3
## 1  2
## 2  1
head(poligonos)
##   id
## 0  3
## 1  2
## 2  1

Se debe notar que por como están programadas estas estructuras espaciales en R, todo supconjunto de un data.frame espacial es nuevamente un data.frame espacial y se pueden generar subconjuntos utilizando los mecanismos usuales para data.frames

# elegimos el primer polígono
poligono <- poligonos[1,]

class(poligono)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(poligono)

Como sucede con un data.frame usual se puede agregar una nueva columna (variable) a un data.frame espacial utilizando el signo “$”.

# agregamos al objetos poligonos la variable intensidad. Este objeto tiene 3 polígonos. Asignarémos un valor de intensidad = 0 al primer polígono, y un valor de intensidad = 2 a los siguientes 2. 
poligonos$intensidad <- c(1,2,2)

head(poligonos)
##   id intensidad
## 0  3          1
## 1  2          2
## 2  1          2
colores <- c("red","blue")

plot(poligonos,col=colores[poligonos$intensidad])

Análisis de capas vectoriales usando R

Existen numerosas librerías para hacer análisis de capas vectoriales de puntos, líneas y polígonos. Se incluye un documento donde se describen algunos de estos paquetes.

Algo que es importante es entender cómo manipular distintas capas vectoriales simultáneamente, por lo que vale la pena leer las funciones disponibles en los paquetes más fundamentales: sp, maptools, raster, etc.

Por ejemplo, es trivial encontrar las intersecciones espaciales entre polígonos y puntos utilizando la función over() del paquete sp.

overlay <- over(puntos,poligonos)

overlay
##    id intensidad
## 1   3          1
## 2   3          1
## 3   3          1
## 4   3          1
## 5   3          1
## 6   3          1
## 7   3          1
## 8   3          1
## 9   3          1
## 10  2          2
## 11  1          2